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Abstract 

The system of non linear reaction-diffusion process in thin 
membrane describing steady state of chemical reactions that 
involves three species is discussed. The equations are 
coupled by the non-linear reaction terms with mixed 
boundary conditions. A closed form of an analytical 
expression of concentrations for the full range of enzyme 
activities has been derived using Homotopy analysis method. 
A simple form of an approximate analytical expression of 
concentrations in terms of dimensionless parameter X is 
also reported. These approximate results are compared with 
the numerical results. A good agreement with simulation 
data is noted. 
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Introduction 

We consider a diffusion controlled chemical reaction 
between two species A and B to form a product, 
according to the reaction mechanism 2A + B — > product . 
This reaction takes place within a thin membrane 
between 'tanks' with abundant supplies of A to the 
left of B to the right of the membrane. We model the 
transport inside the membrane as diffusive, thus the 
model will be given by a system of non-linear 
reaction-diffusion equations that are coupled with the 
non-linear reaction terms. The reaction path consists of 
a coupled pair of rapid irreversible simple reaction 
mechanism [Ariel (2010)]. 



A + B- 



->C, 



A + C — — > product 



(1) 



Seidman et. al (1997 and 2005) and Kalacheve et. al. 
(2003) provided the rigorous singular perturbation 
analysis for the steady-state problem. The 
corresponding non-steady state system of this problem 
has been considered by Haario Seidman (1994), for the 
complex boundary conditions, to describe reactions in 
the film model for a gas/liquid interface. Also the 
steady state problem has many important applications, 
in chemical engineering modeling. Recently, Butuzov 
et al. (1999 and 2007) have obtained some related 
problems of exchange of stabilities using different 
techniques (upper and lower solutions). However, to 
the best of author's knowledge, no general analytical 
results of substrate concentration for all values of 
dimensionless parameter X have been published. The 
purpose of this communication is to derive 
approximate analytical expressions for the steady-state 
concentrations for all values of X using Homotopy 
analysis method. 

Mathematical Formulation of the Problems 

The governing non-linear reaction diffusion equation 
in a thin membrane is expressed in the following non- 
dimensional format [Ariel (2010)]: 



u xx = Xuv + uw 
Vxx = Xuv 
w„ =uw- Xuv 



(2) 
(3) 
(4) 



where u(x), v(x) and w(x) denote the concentrations of 
the chemical species A, B and C respectively. The 
diffusion coefficients of three species are considered to 
have an equal diffusion coefficient which is equal to 1. 
We assume that the specie A is supplied with a given 
fixed concentration a > at x = 0, and the specie B 
with j3 > at x = 1. Boundary conditions are 



where X and ju denote the binary reaction rates. 



-a; v x = ; w = y at x =0 



(5) 



10 



International Journal of Automation and Control Engineering Volume 2 Issue 1, February 2013 



www.seipub .org/ij ace 



Q;v = /3; w x =0 at x = 1 



(6) 



Since the appearance of the large factor X is much 
greater than 1 in one of the terms in each reaction- 
diffusion equations (2)-(4) the equations have 
singularly perturbed problems [Ariel (2010)]. 

Solution of Boundary Value Problem Using 
Homotopy Analysis Method (HAM) 

Liao (1992, 1995, 2003, 2004, 2010 and 2012) proposed 
a powerful analytical method for non-linear problems, 
namely the Homotopy analysis method. This method 
provides an analytical solution in terms of an infinite 
power series. However, there is a practical need to 
evaluate this solution and to obtain numerical values 
from the infinite power series. In order to investigate 
the accuracy of the Homotopy analysis method (HAM) 
solution with a finite number of terms, the differential 
equations in the system were solved. The Homotopy 
analysis method [Liao (1992, 1995, 2003, 2004, 2010, 
2012), Eswari et.al (2010) and Jafari et. al (2009)] is a 
preferred technique comparing to another 
perturbation method. 

Homotopy perturbation method [Chowdhury et. 
al (2007), Eswari et.al (2010), Ghori et. al (2007), Coyle 
et. al (1986), Ozis et. al (2007), Madden et. al (2003) ] is 
a special case of Homotopy analysis method. Different 
from all reported perturbation and non perturbative 
techniques, the Homotopy analysis method itself 
provides us with a convenient way to control and 
adjust the convergence region and rate of 
approximation series, when necessary. Briefly 
speaking, the Homotopy analysis method has the 
following advantages: itt is valid even if a given 
nonlinear problem does not contain any small/large 
parameter at all; it can be employed to efficiently 
approximate a nonlinear problem by choosing 
different sets of base functions. The Homotopy 
analysis method contains the auxiliary parameter h 
which provides us with a simple way to adjust and 
control the convergence region of solution series. 
Using this method, we can obtain the following 
solution to (2) - (6) (see Appendix B). 



w(x) = ha(y - X/3)x + 



u(x) = a + ha\X/3 + y] 



v(x) = /3 + 



and 



hXaj3(l-x 2 )^ 



(7) 



(8) 



ha{X/3 - ay)x z 



(9) 



The equations (7)-(9) represent the new analytical 
expression of concentration of species for all values of 
dimensionless parameter. The reaction rate q is given 
by 



q = Xuv = X 



a + ha[ip + y] 



x^ 



hXaj3(\-x 1 )' 



(10) 



Numerical Simulation 

In order to find reveal the accuracy of our analytical 
method, the non-linear differential eqns. (2)-(4) are 
also solved by numerical methods. The function bup4c 
in Matlab software which is a function of solving two- 
point boundary value problems (TPBVPs) for ordinary 
differential equations is used to solve these equations 
numerically. Our analytical result is compared with 
numerical solution and it gives a satisfactory 
agreement (See figures (1-3)). The Matlab program is 
also given in Appendix D. 

Results and Discussions 

Fig. 1 represents the normalized steady-state 
concentration u(x) for different values of 
dimensionless parameter X = 0.1, 0.5, 1, 3 and 5. From 
this figure, it is evident that the values of the 
concentration decreases when dimensionless 
parameter X increases to x = l. Fig. 2 shows the 
normalized steady-state concentration v(x) versus the 
dimensionless distance x for various values of 
dimensionless parameter X . From this figure, it is 
obvious that the values of the concentration increases 
when dimensionless parameter A decreases to x = . 
The normalized steady-state concentration w(x) 
versus the dimensionless distance x for various values 
of dimensionless parameter X is plotted in the Fig. 3. 
From this figure, it is inferred that the value of the 
concentration will increase, when the diffusion 
parameter X increases. Fig.4 shows the normalized 
steady state concentrations u(x), v(x) and w(x) for some 
fixed value of X = 5 . From this figure it is observed 
that the concentration u(x) of species A from its initial 
value, where the concentrations v(x) of the species B 
and w(x) of the species C increase from its initial 
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values. 




0.4 I 1 1 1 1 1 1 1 1 1 1 

0.1 2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Dimensioiiless distance x 

Fig. 1 Normalized steady-state concentration u{x) versus 
the dimensionless distance x. the concentrations were 
computed using (7) for various values of the dimensionless 

PARAMETER A AND h = -0.2. 



Fig. 3 Normalized steady-state concentration w(x) versus 
the dimensionless distance x . the concentrations were 
computed using (9) for various values of the dimensionless 

PARAMETER A AND h = -0.2. 




Dimensioiiless distance x 




. ■ ■ ' Analytical 

*° *********** Numerical 



cr=1.6,/3=O.B, y = 0.01 



0.1 0.2 0.3 0.4 OJ M 0.7 0.8 0.9 1 
Dimensioiiless distance X 



Fig. 4. Normalized steady-state concentration 

U(x), \'(x) AND w(x) VERSUS THE DIMENSIONLESS DISTANC X . 
THE CONCENTRATIONS WERE COMPUTED USING (7) - (9) AND FOR THE 
FIXED VALUE OF THE DIMENSIONLESS PARAMETER A - 5 AND 

h = -0.2. 

Figs. 5-7 shows the dimensionless reaction rate q 
using for various values of A . Thus, it is concluded 
that there is a simultaneous increase in the values of 
the reaction rate as well as in A for the fixed value of 
a,p and y . The optimum value of h can be obtained 
using h curve which is given in the Figures 8-9. 



Fig. 2 Normalized steady-state concentration v(x) versus the 
dimensionless distance x . the concentrations were computed 
using (8) for various values of the dimensionless parameter x 
AND h = -0.2. 





0.1 0.2 0.3 0.4 0.5 0,6 0.7 0.8 0.9 1 
Dimensionless distance x 

FIG. 5 DIMENSIONLESS REACTION RATE VERSUS THE DIMENSIONLESS 
DISTANCE X USING (10) FOR THE VALUE OF THE DIMENSIONLESS 
PARAMETER A = 0.1 WHEN a = 1.6, ft - 0.8, / - 0.01 AND 

h = -0.2. 



12 



International Journal of Automation and Control Engineering Volume 2 Issue \, February 2013 



www.seipub.org/ijace 



1.2 



1 18 - 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
D im eitsio nless distance x 



FIG. 6 DlMENSIONLESS REACTION RATE VERSUS THE DIMENSIONLESS 
DISTANCE X USING (10) FOR THE VALUE OF THE DIMENSIONLESS 
PARAMETER X = 1, WHEN a = 1.6, /3 = 0.8, ^" = 0.01 AND 

/1 = -0.2. 



4 




1 1 1 ' ' 1 ' 1 ' 1 1 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 8 0.9 1 

D im ensio nless distances 

FIG. 7 DIMENSIONLESS REACTION RATE VERSUS THE DIMENSIONLESS 
DISTANCE X USING (10) FOR THE VALUE OF THE DIMENSIONLESS 
PARAMETER 1=5, WHEN a =1.6, /? = 0.8, ^" = 0.01 AND 

h = -0.2. 




-1.2 - 



-1,4 ' ' ' 1 1 ' 1 ' 1 

-0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 

h 

Fig. 8 The h curve to indicate the convergence region u (0. 1) 
for when a = 1 .6, /? = 0.8, 7 = 0.01 and A -5. 




1.46- 

-0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 
h 

Fig. 9 The h curve to indicate the convergence region m(0. 1) 
for when a = 1 .6, /? = 0.8, / = 0.01 and X =5. 

Conclusion 

The system of time independent reaction-diffusion 
equation coupled through the non linear reaction 
terms in thin membrane has been solved analytically 
and numerically. Analytical expressions of the 
concentrations of species are derived by using the 
Homotopy analysis method. The primary result of this 
work is simple and approximate expressions of the 
concentrations for all values of the dimensionless 
parameter X . This analytical result will be useful to 
analyze the behavior of the internal layers. This 
method is an extremely simple and it is also a 
promising method to solve other non-linear equations. 
This method can be easily extended to find the 
solution of all other non-linear equations. 
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Appendix A: Basic Concept of Homotopy Analysis 
Method (HAM) 

Consider the following differential equation: 

M«(f)] = (A.l) 

Where N is a nonlinear operator, t denotes an 
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independent variable, u(t) is an unknown function. For 
simplicity, we ignore all boundary or initial conditions, 
which can be treated in the similar way. By means of 
generalizing the conventional Homotopy method, 
Liao (2012) constructed the so-called zero-order 
deformation equation as: 



(1 - p)Utfr, p) - u Q (f )] = P hH(t)N[cp(t; p)] 



(A.2) 



where pe [0,1] is the embedding parameter, h * is a 
nonzero auxiliary parameter, H(t) * is an auxiliary 

function, L an auxiliary linear operator, U ° (t) is an 
initial guess of u(t), <p(t : p) is an unknown function. It 
is important to note that one has great freedom to 
choose auxiliary unknowns in HAM. Obviously, when 
p = and p = 1 , it holds: 

cp{ffi)=u {t) and cp{fX) = u{t) (A.3) 
respectively. Thus, as p increases from to 1, the 
solution (p(t; p) varies from the initial guess u (t) to 
the solution u (f). Expanding <p(t; p) in Taylor series 
with respect to p, we have: 



<p(t;p) = u Q (t) + Y J u m (t)p" 
where 



m=l 



1 d*«t,ph 



(A.4) 



(A.5) 



If the auxiliary linear operator, the initial guess, the 
auxiliary parameter h, and the auxiliary function are 
so properly chosen, the series (A.4) converges at p =1 
then we have: 



u(t) = u Q (t) + ^u m (t) 



(A.6) 



Applying L 1 on both side of equation (A7), we get 



«m (0 = Zm*m-\ (0 + hL- 1 [ff(f )M n {U m _ x )] 



(A10) 



In this way, it is easily to obtain u m for m>\, at M th 
order, we have 



M 

«(0 = ^K m (0 



(A.ll) 



m=0 



WhenM -»+oo, we get an accurate approximation of 
the original equation (A.l). For the convergence of the 
above method we refer the reader to Liao [20]. If 
equation (A.l) admits unique solution, then this 
method will produce the unique solution. 

Appendix B: Solution of Non-linear Equations (2) to (6) 
Using HAM 

In this Appendix, we indicate how (7) to (9) in this 
paper are derived. To find the solution of (2), (3) and 
(4), it can be simplified to [Ariel (2010)] 



j2 

a u 

dx 2 



■ A-uv - uw = 



, 2 

a v 
dx 2 



Auv =0 



d 2 w 
dx 2 



-uw + Xuv = 



Differentiating (A.2) for m times with respect to the 
embedding parameter p, and then setting p = and 


(1- 


P) 


d 2 u 
dx 2 _ 


= hp 


d 2 u 
dx 2 


-Auv 


-uw 


finally dividing them by ml, we will have the so-called 
















mth -order deformation equation as: 


(1- 


P) 


d 2 v 
dx 2 _ 


= hp 


d v 

dx 2 


-Auv 




Uu m -X m u m -x\ = hH(m m (u m -i) (A.7) 
















where 


(1- 


P) 


d 2 w 
dx 2 


= hp 


d 2 w 
dx 2 


+ Auv 


-uw 



9t_(K m -l) = 



1 d m - l N[p(r,p)] 

(m-iy. d P m - 1 



(B.l) 
(B.2) 

(B.3) 

(B.4) 
(B.5) 

(B.6) 
(B.7) 
(B.8) 

(A.8) -j^g a pp rox i ma te solution of (B.l) and (B.2) and (B.3) is, 



Now the boundary conditions becomes 

x=0, u = a , — = 0, w = y 
dx 

du dw 
x=l, — = Q,v = /3 , — = 

dx dx 

We construct the Homotopy as follows 



And 

X m ~ 



U = M Q + pU x + p U 2 +. 



[0, m<l, 
ll, m>\. 



(A.9) 



V = V Q + pVy + p V 2 +- 



w = w Q + pw x + p w 2 +- 



(B.9) 

(B.10) 
(B.ll) 
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The initial approximations are as follows 




u (0) = a and w (l) = 


(B.12) 


m ; (0) = and w ( '(l) = , i = 1,2... 


(B.13) 


v (0) = and v (l) = ^ 


(B.14) 


v;(0) = and v,(l)=0, i = 1,2... 


(B.15) 


w (0) = 7 and w (l) = 


(B.16) 


w t (0) = and w] (1) = 0, i = 1,2... 


(B.17) 


Substituting (B.9) to (B.ll) into (B.6) to (B.8) we have 



(1-p) 



d 2 (u + pu t +...) 



dr 1 



= hp 



d (u + pii\ +...) 

dx 2 

-Muo+pii! + )(v +pv { + ) 

-(uo+pui +....)(w +pw t +....) 



(B.18) 



a-p) 



d (v +pv x + ..) 



dx 2 



= hp 



d (v +pv x + ..) 

dx 2 

-A(u +pu 1 + ....)(v +pv l + ) 



(B.19) 



d-p) 



d {wq + pw\ + ....) 

dx 1 



-hp 



d {wq + pw\ + ....) 

dx 2 

+ Muq + pu\ + ...)(vq + pv\ + ...) 

- (UQ + pU\ + ....)(\VQ + p\V\ + ...). 



(B.20) 



Comparing the coefficients of like powers of p in (B.18) 
we get 

o d Ua 



dx 2 



= 



P ■ ^-(h + l) 

dx 



d 2 u a 
dx 1 



(B.21) 

+ hAu v + hu w = (B.22) 



Comparing the coefficients of like powers of p in (B.19) 
we obtain. 



,2 

o.d v 



dx 1 



= 



(B.23) 



dx 2 



(h + 1) 



,2 

d v 
dx 2 



+ hA.u Q v Q =0 



(B.24) 



Comparing the coefficients of like powers of p in (B.20) 
we have 



= 



dx' 



^2 
\ d u\ 

P : ^--(h + l) 

dx 1 



d WQ 
dx 2 



(B.25) 

+ huQWQ - h/luQVQ = (B.26) 



Solving (B.21) to (B.26) and using the boundary 
conditions (B.12) to (B.17), we can obtain the following 
results: 



uq =a, u\ = ha\xp + y\x , vq = (3, 

n 



' hXaPil-x 1 )^ 



wq =0, w\= ha(y - Xf5)x + 



ha(Xp - ay)x 



According to the HPM, we can conclude that 
u = lim u(x) = u Q + Wj 
v = lim v(x) = v + Vj 

w = lim w(x) = w + w l 
P ->\ 



(B.27) 
(B.28) 

(B.29) 
(B.30) 

(B.31) 



After putting (B.27) and (B.28) into (B.29) to (B.31), we 
obtain (7) - (9) in the text. 

Appendix C: Determining the Region of h for Validity 

The analytical solution should converge. It should be 
noted that the auxiliary parameter h controls the 
convergence and accuracy of the solution series. The 
analytical solution represented by (7) contains the 
auxiliary parameter h which gives the convergence 
region and rate of approximation for the Homotopy 
analysis method. In order to define region such that 
the solution series is independent of h, a multiple of 
h curves are plotted. The region where the 
concentration u(x) and u\x) versus h is a horizontal 

line known as the convergence region for the 
corresponding function. The common region among 
u(x) and its derivatives are known as the over all 
convergence region. To study the influence of h on 
the convergence of solution, h curves of w(0.1) and 
m'(0.1) are plotted in Fig. (7) and (8) respectively for 
a = 1.6, ^ = 0.8,^ = 0.01 and A = 5. These figures 
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clearly indicate that the valid region of h is about (-0.4 
to -0.05). Similarly we can find the value of the 
convergence control parameter h for different values 
of constant parameters. 

Appendix D : Matlab Program to Find the Numerical 
Solution of Non-linear Equations (2) to (6): 

function pdex4 

m = 0; 

x = linspace(0,l); 
t=linspace(0, 100000); 

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,X/t); 

ul = sol(:,:,l); 

u2 = sol(:,:,2); 

u3=sol(:,:,3); 

figure 

plot(x,ul(end,:)) 
title('ul(x,t)') 
xlabel( 'Distance x') 
ylabel('ul(x,2)') 

% 

Figure 

plot(x,u2(end,:)) 
title('u2(x,t)') 
xlabel('Distance x') 
ylabel('u2(x,2)') 



Figure 

plot(x,u3(end,:)) 
title('Solution at t = 2') 
xlabel( 'Distance x') 

ylabel('u3(x,2)') 

% 



function [c,f,s] = pdex4pde(x,t,u,DuDx) 

c = [l; 1;1]; 

f=[l; 1;1] *DuDx; 

y = u(l)*u(2); 

yl=u(l)*u(3); 



alpha=1.6; 

gamma=0.01; 

beta=0.8; 

lamta=5; % parameters 
F =(-lamta*y-yl); 

Fl=(-lamta*y); % non linear terms 

F2=(lamta*y-yl); 

s=[F;Fl;F2]; 

% 

Function uO = pdex4ic(x); 
%create a initial conditions 
uO = [0; 1; 0]; 

/ G 

function[pl,ql,pr,qr]=pdex4bc(xl / ul,xr,ur,t) 
%create a boundary conditions 

pl=[ul(l)-1.6; 0;ul(3)-0.01]; 

ql = [0; 1; 0]; 

pr = [0;ur(2)-0.8;0]; 

qr = [l;0;l]; 

Appendix E : Nomenclature 



Symbol 


Meaning 


U 


Concentration of the chemical species A 


V 


Concentration of the chemical species B 


W 


Concentration of the chemical species C 


X 


Dimensionless parameter 


X 


Dimensionless distance 


a 


Fixed concentration of the species A 


P 


Fixed concentration of the species B 


7 


Fixed concentration of the species C 


q 


Dimensionless reaction rate 
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